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ABSTRACT 

We study the interplay of clumping at small scales with the collapse and relaxation of per- 
turbations at much larger scales. We present results of our analysis when the large scale per- 
turbation is modelled as a plane wave. We find that in absence of substructure, collapse leads 
to formation of a pancake with multi-stream regions. Dynamical relaxation of plane wave is 
faster in presence of substructure. Scattering of substructures and the resulting enhancement 
of transverse motions of haloes in the multi-stream region lead to a thinner pancake. In turn, 
collapse of the plane wave leads to formation of more massive collapsed haloes as compared 
to the collapse of substructure in absence of the plane wave. The formation of more massive 
haloes happens without any increase in the total mass in collapsed haloes. A comparison with 
the Burgers' equation approach in absence of any substructure suggests that the preferred 
value of effective viscosity depends primarily on the number of streams in a region. 
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1 INTRODUCTION 

Large scale structures like galaxies and clusters of galaxies are 
believed to have formed by gravitational amplification of small 
perturbations (Peebles, 1980; Peacock, 1998; Padmanabhan, 2002; 
Bernardeau et al., 2002). Observations suggest that the initial den- 
sity perturbations were present at all scales that have been probed 
by observations. An essential part of the study of formation of 
galaxies and other large scale structures is thus the evolution of 
density perturbations for such initial conditions. Once the ampli- 
tude of perturbations at any scale becomes large, i.e., 5 ~ 1, the 
perturbation becomes non-linear and the coupling with perturba- 
tions at other scales cannot be ignored. Indeed, understanding the 
interplay of density perturbations at different scales is essential for 
developing a full understanding of gravitational collapse in an ex- 
panding universe. The basic equations for this have been known for 
a long time (Peebles, 1974) but apart from some special cases, few 
solutions are known. 

A statistical approach to this problem based on pair conserva- 
tion equation has yielded interesting results (Hamilton et al., 1991; 
Nityananda and Padmanabhan, 1994; Padmanabhan, 1996; Engi- 
neer, Kanekar and Padmanabhan, 2000), and these results have 
motivated detailed studies to obtain fitting functions to express the 
non-linear correlation function or power spectrum in terms of the 
linearly evolved correlation function (Hamilton et al., 1991; Iain, 
Mo and White, 1995; Peacock and Dodds, 1996; Smith et al., 
2003). 

It is well known from simulation studies that at the level of 
second moment, i.e., power spectrum, correlation function, etc.. 



large scales have an important effect on small scales but small 
scales do not have a significant effect on large scales (Peebles, 
1985; Klypin and Melott, 1992; Little, Weinberg and Park, 1991; 
Bagla and Padmanabhan, 1997a; Couchman and Peebles, 1998). 
Most of these studies used power spectrum as the measure of clus- 
tering. Results of these simulation studies form the basis for the use 
of N-Body simulations, e.g., from the above results we can safely 
assume that small scales not resolved in simulations do not effect 
power spectrum at large scales and can be ignored. 

Substructure can play an important role in the relaxation pro- 
cess. It can induce mixing in phase space (Lynden-Bell, 1967; 
Weinberg, 2001), or change halo profiles by introducing transverse 
motions (Peebles, 1990; Subramanian, 2000), and, gravitational in- 
teractions between small clumps can bring in an effective collision- 
ality even for a coUisionless fluid (Ma and Boylan-Kolchin, 2004; 
Ma and Bertschinger, 2003). Thus it is important to understand the 
role played by substructure in gravitational collapse and relaxation 
in the context of an expanding background. In particular, we would 
like to know if this leaves an imprint on the non-linear evolution of 
correlation function. Effect of substructure on collapse and relax- 
ation of larger scales is another manifestation of mode coupling. 

In this paper, we report results from a study of mode coupling 
in gravitational collapse. In particular, we study how presence of 
density perturbations at small scales influences collapse and relax- 
ation of perturbations at larger scales. These effects have been stud- 
ied in past (Evrard and Crone, 1992) but the motivation was slightly 
different (Peebles, 1990). We believe it is important to study the is- 
sue in greater detail and make the relevance of these effects more 
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quantitative using N-Body simulations with a larger number of par- 
ticles. We also study the reverse process, i.e., how does collapse of 
perturbations at large scales effect density perturbations at much 
smaller scales. 

It is well known that the local geometry of collapse at the time 
of initial shell crossing is planar in nature (Zel'dovich, 1970), hence 
we model density perturbations as a single plane wave in this work. 
Simple nature of the large scale fluctuation allows us to study inter- 
action of well separated scales without resorting to statistical esti- 
mators like power spectrum. We are studying the same problem in 
a more general setting and those results will be reported in a later 
publication. 

Key features of collapse of a plane wave can be understood us- 
ing quasi-linear approximations, at least at a qualitative level. Initial 
collapsing phase is well modelled by the Zel'dovich approxima- 
tion (Zel'dovich, 1970), wherein particles fall in towards the centre 
of the potential well. Zel'dovich approximation breaks down af- 
ter orbit crossing as it does not predict any change in the direction 
of motion for particles, thus in this approximation particles con- 
tinue to move in the same direction and the size of the collapsed 
region grows monotonically. In a realistic situation we expect par- 
ticles to fall back towards the potential well and oscillate about 
it with a decreasing amplitude, and the collapsed region remains 
fairly compact. Several approximations have been suggested to im- 
prove upon the Zel'dovich approximation (Gurbatov, Saichcv and 
Shandarin, 1989; Shandarin and Zeldovich, 1989; Weinberg and 
Gunn, 1990; Matarrese et al., 1992; Brainerd, Scherrer and Villum- 
sen, 1993; Bagla and Padmanabhan, 1994; Sahni and Coles, 1995; 
Hui and Bertschinger, 1996). The adhesion approximation (Gur- 
batov, Saichev and Shandarin, 1989; Weinberg and Gunn, 1990) 
invokes an effective viscosity: this prevents orbit crossing and con- 
serves momentum to ensure that pancakes remain thin and matter 
ends up in the correct region. This changes the character of motions 
in dense regions (no orbit crossing or mixing in the phase space) but 
predicts locations of these regions correctly. If one assumes that the 
gravitational potential evolves at the linear rate (Brainerd, Scherrer 
and Villumsen, 1993; Bagla and Padmanabhan, 1994), then it can 
be shown that the collapsed region remains confined. The effective 
drag due to expanding background slows down particles and they 
do not have enough energy to climb out of the potential well. 

Thus the process of confining particles to a compact collapsed 
region results from a combination of expansion of the universe and 
gravitational interaction of in falling particles. None of the approx- 
imations captures all the relevant effects. Therefore we must turn to 
N-Body simulations (Bertschinger, 1998; Bagla, 2004) in order to 
study collapse and relaxation of perturbations in a complete man- 
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2 EVOLUTION OF PERTURBATIONS 

We will consider only gravitational effects here and ignore all other 
processes. We assume that the system can be described in the New- 
tonian limit. The growth of perturbations is then described by the 

coupled system of Euler's equation and Poisson equation in comov- 
ing coordinates along with mass conservation, e.g., see (Peebles, 
1980). 
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It is assumed that the density field is generated by a distribution 
of particles, each of mass rrii, position Xj. Ho is the present value 
of Hubble constant, Qnr is the present density parameter for non- 
relativistic matter and a is the scale factor. In this paper wc will 
consider the Einstein-de Sitter universe as the background, i.e., 
fl„r = 1. These can be reduced to a single non-linear differential 
equation for density contrast (Peebles, 1974). 
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The terms A and B are the non-linear coupling terms between dif- 
ferent modes. B couples density contrasts in an indirect manner 
through velocities of particles (x,). The equation of motion still 
needs to be solved for a complete solution of this equation, or we 
can use some ansatz for velocities to make this an independent 
equation. 

It can be shown that individual virialised objects, i.e., objects 
that satisfy the condition 2T+U = where T is the kinetic energy 
and U is the potential energy, do not make any contribution towards 
growth of perturbations through mode coupling (Peebles, 1974) at 
much larger scales, i.e., the A — B term is zero. The contribution 
of mode coupling due to interaction of such objects is not known. 

Approximate approaches to structure formation can be devel- 
oped by ignoring interaction of well separated scales. The evolu- 
tion of density perturbations can be modelled as a combination of 
non-linear collapse at small scales, and the collapsed objects can 
be displaced using quasi-linear approximations (Bond and Myers, 
1996a; Bond and Myers, 1996b; Bond and Myers, 1996c; Monaco 
et al., 2002; Monaco, Theuns and Taffoni, 2002; Taffoni, Monaco 
and Theuns, 2002). These approaches yield an acceptable descrip- 
tion of properties of collapsed objects and their distribution for a 
first estimate. PINOCCHIO (Monaco et al., 2002; Monaco, Theuns 
and Taffoni, 2002; Taffoni, Monaco and Theuns, 2002) provides 
sufficient information about halo properties and merger trees for 
use with semi-analytic models of galaxy formation. The efficacy of 
these models puts an upper bound on the effects of mode coupling 
that we are studying here. 

In this paper we simplify the system by starting with pertur- 
bations that have a non-zero amplitude only for two sets of scales. 
We simulate the collapse of a plane wave by starting with non-zero 
amplitude of perturbations for the fundamental mode of the sim- 
ulation box along the z axis, the wavenumber of the fundamental 
mode is denoted by fc/ . This serves as the large scale perturbation 
in our study. The amplitude for this mode is chosen so that shell 
crossing takes place when the scale factor a — 1. Power spectrum 
for small scale fluctuations is chosen to be non-zero in a range of 
wave numbers ko ± Sk with a constant amplitude across this win- 
dow, i.e., As(fc) — a A for ko — Sk ^ k ^ ko + Sk. A Gaussian 
random realisation of this power spectrum is used for small scale 
fluctuations. Here A'l{k) is the power per logarithmic interval in 
k contributed by small scales (large k) and A is the amplitude of 
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Table 1. This table lists parameters of N-Body simulations we have used. 
All the simulations used 128^ particles. The first column lists name of the 
simulation, second column Usts the code that was used for running the sim- 
ulation, third column gives the relative amplitude of small scale power and 
the plane wave, the fourth column tells us whether the large scale plane 
wave was present in the simulation or not, and the last column lists the dis- 
tribution of particles before these are displaced using a realisation of the 
power spectrum. Grid distribution means that particles started from grid 
points. Perturbed grid refers to a distribution where particles are randomly 
displaced from the grid points, this displacement has a maximum ampUtude 
of 0.05 grid points. Such an initial condition is needed to prevent particles 
from reaching the same position in plane wave collapse as such a situation is 
pathological for tree codes. The TreePM simulations were run with a force 
softening length equal to the grid length. 
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Figure 1. Phase space plot for the plane wave at late stages of collapse 
for the simulation PM.OOL. Velocity of particles is plotted as a function of 
position. Regions where particles with different velocities can be found are 
the multi-stream regions. As we approach the centre of pancake located at 
z = 32, we go from a single stream region to a three stream region and so 
on up to seven stream region near the centre. 
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the fundamental mode that gives rise to the plane wave. The ratio 
of A'^ (k) at k — and for the plane wave is denoted by q, thus 
when a = 1 collapse of perturbations at these scales happens at the 
same time whereas for a > 1 perturbations at small scales collapse 
before the plane wave collapses. We chose the ratio fco/fc/ = 8 so 
that there is distinct separation in the scales involved. 

Collapse of a plane wave of this kind leads to formation of a 
multi-stream region, we will also use the term pancake to describe 
this region. Figure 1 shows the phase space plot for the plane wave 
at late stages of collapse for the simulation PM.OOL (see table 1 
for details of the simulation). Velocity of particles is plotted as a 
function of position, only the z component is plotted as there is no 
displacement or velocity along other directions in this simulation. 
Regions where particles with different velocities can be found are 
the multi-stream regions. As we approach the centre of pancake 
located at z = ,32, we go from a single stream region to a three 
stream region and so on up to seven stream region near the centre. 

In initial stages, the mass in the pancake increases rapidly as 
more particles fall in. Figure 2 shows this in terms of over-density 
which increases sharply from a = ltoa = 2.A significant fraction 
of the total mass falls into the pancake and the infall velocities for 
the remaining matter are very small. In this regime the mass of 
the pancake is almost constant, this can be seen from the panels of 
figure 8 where the mass enclosed in the pancake region is almost 
constant from a = 2 to a = 4. 

In absence of any substructure the coUisionless collapse re- 
tains planar symmetry and we have layers of multi-stream regions 
with the number of streams increasing towards the centre of the 
pancake. Presence of small scale fluctuations can induce transverse 
motions and these motions are amplified in the pancake. 

Weakly bound substructure can be torn apart due to interaction 
with rapidly in falling matter. On the other hand, higher average 
density in the multi stream region can lead to rapid growth of per- 



turbations. It is known that pancakes are unstable to fragmentation 
due to growth of perturbations (Valnia et al., 1997). The velocity 
field is anisotropic due to infall along one direction, hence the rate 
at which perturbations grows will also exhibit anisotropics. Veloc- 
ity dispersion along the direction of plane wave collapse is larger 
than the transverse direction, hence the growth of fluctuations in 
the transverse plane is expected to be more rapid. 

If the in falling material contains collapsed substructure, then 
gravitational interactions between these can induce large transverse 
velocities. This takes away kinetic energy from the direction of in- 
fall, which in turn can lead to more fragmented and thinner multi- 
stream region. 

In the following sections we describe the numerical experi- 
ments have undertaken in detail, and test the physical ideas and 
expectations outlined above. 



3 NUMERICAL EXPERIMENTS AND RESULTS 

We used a Particle-Mesh code (Bagla and Padmanabhan, 1997b) 
and the TreePM code (Bagla, 2002; Bagla and Ray, 2003). Some 
simulations were run using the parallel TreePM (Ray and Bagla, 
2004). TreePM simulations used spline softening with softening 
length equal to the length of a grid cell in order to ensure coUi- 
sionless evolution. We used force softening assuming a spline ker- 
nel (Springel, Yoshida and White, 2001). All the simulations were 
carried out with 128'' particles. Table 1 lists parameters of the sim- 
ulations we have used for this paper. We have used two types of 
initial distribution of particles. In the Grid distribution particles are 
located at grid points before being displaced to set up the initial per- 
turbations. Perturbed grid refers to a distribution where particles are 
randomly displaced from the grid points (Bagla and Padmanabhan, 
1997b), this displacement has a maximum amplitude of 0.05 grid 
length. Such an initial condition is needed to prevent particles from 
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Figure 2. Density profile is plotted for two epochs for simulations PM.OOL 
and T.OOL - in these simulations density varies only in the direction along 
the plane wave. SoUd lines show the density profile at a = 1 and a = 2 
from the PM.OOL simulation, dashed lines show the density profile from the 
T.OOL simulation with the same profile. 

reaching the same position in plane wave collapse as such a situ- 
ation is pathological for tree codes. These small displacements do 
not affect the power spectrum to be realised, PM.OOL and T.OOL 
were compared to test for any systematic effects. 

Simulations T.lOP and T.40P were similar to T.IOL and 
T.40L except that the small scale fluctuations were restricted to 
the direction orthogonal to the direction of plane wave. Thus the 
small scale fluctuations had the same form for all z. These simula- 
tions are useful for differentiating between competing explanations 
for results outlined below. 

In addition to the N-Body simulations listed in table 1, we also 
carried out one dimensional simulations within the adhesion model 
(Gurbatov, Saichev and Shandarin, 1989) with a finite viscosity fol- 
lowing a method similar to the one outlined by Weinberg and Gunn 
(1990). 

Figure 2 shows the density profile of the pancake for PM.00L 
and T.OOL simulations at two epochs. These figures demonstrate 
that the density profiles in these simulations are almost identical, 
indeed the tiny differences can be attributed to the different initial 
distribution of particles. We have checked this assertion by miming 
the PM.OOL with the perturbed grid initial conditions. The TreePM 
method has a slightly better resolution but it does not induce any 
new features. This is expected as the force softening length used 
in the TreePM simulations is one grid length, same as the aver- 
age inter-particle separation and it has been shown than such force 
softening does not induce two body collisions (Splinter et al., 1998; 
Melott et al., 1997). We will mostly use TreePM simulations for the 
remaining part of this study. 



3.1 Thickness of pancake 

An important indicator of the role played by substructure is the 
thickness of the pancake that forms by collapse of the plane wave. 
If substructure does not play an important role in evolution of large 
scale perturbations then the thickness of pancake should not change 
by a significant amount. On the other hand, if substructure does 
indeed speed up the process of dynamical relaxation then we should 
see some signature in terms of the thickness of pancake, velocity 
structure, or both. Any such effect will be apparent only at late 
times as infall of matter into the pancake dominates at early times. 



Dynamical effects of substructure will become important only at 
late times. 

Figure 3 shows a slice from some of the simulations listed in 
table 1. The plane wave coUapses along the vertical axis. Config- 
uration at a = 2 is shown here, the plane wave begins to collapse 
at a = 1. Different panels in this figure refer to simulation T.OOL, 
T.IOL and T.40L. The boundary of the multi stream region is vis- 
ible clearly in all the slices even though this region is fragmented 
in the last panel (T.40L). It is clear that the pancake is thinner in 
simulations with more substructure. 

A more detailed comparison of simulations with different 
level of substructure is shown in figure 4. The top panel of this 
figure shows the averaged over-density as a function of the z co- 
ordinate, the plane wave collapses along this axis. Over-density is 
averaged over all x and y for a given interval (z ± Az) to obtain 
the averaged values plotted here. The peak over-density at the cen- 
tre of the pancake is smaller in simulations with more substructure. 
The mass enclosed within a given distance of the centre of pancake 
(defined here as the trough of the potential well of the plane wave) 
is smaller for more substructure, even though the variation is very 
small at less than 10% between the extreme cases (see figure 8). 
Potential wells corresponding to substructure prevent infall into the 
pancake region. As the amount of substructure is increased, there is 
visible reduction in the size of the region around the pancake where 
density is greater than average. The visual impression of figure 3 is 
reinforced by the variation of over-density. 

The middle panel of figure 4 shows the rms velocities of parti- 
cles in direction transverse to the plane wave collapse as a function 
of the z coordinate. As in the top panel, averaging is done over all 
X and y for a given interval {z ± Az). This plot shows that the 
transverse motions are enhanced in the dense pancake region. The 
amplitude of transverse motions is larger in simulations with more 
substructure. Size of the region where these motions are significant 
varies with the amount of substructure, as in case of over-density 
(top panel). The rms transverse velocities do not go to zero outside 
the pancake region, instead these level off to a small residual value. 

Transverse motions are due to motions of particles in clumps 
that constitute substructure, due to infall of particles in these 
clumps, and, transverse motions of clumps as they move towards 
each other. In order to delineate these effects, we have plotted 
the rms velocities for haloes in the last panel of figure 4. These 
haloes were selected with the friends-of-friends (FOF) algorithm 
using a Unking length of Z = 0.2 grid length. Transverse compo- 
nent of the velocity of centre of mass for haloes with more than 
50 particles was used for the figure. Such a high cutoff for halo 
members is acceptable because typical haloes have several hundred 
members, see the following subsection on mass functions. Differ- 
ences between simulations with different amount of substructure 
are more pronounced than in the middle panel. For simulations with 
a small amount of substructure, motion of clumps is subdominant 
and hence the transverse motions are contributed mostly by internal 
motions and infall. In simulations with more substructure, motions 
of clvunps contribute significantly to the rms transverse velocity. 
Gravitational attraction of clumps, particularly in close encounters 
in the pancake region induce the transverse component. Collisions 
are enhanced in the pancake region as the number density of clumps 
is higher. 

In order to convince ourselves that transverse motions induced 
by scattering/collision of clumps is the most likely reason for the 
reduced thickness of pancakes, we compare simulations T.IOL and 
T.40L with T.lOP and T.40P. In T.lOP and T.40P simulations, 
the small scale fluctuations do not have any z dependence. In these 
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Figure 3. Panels in this figure sliow the same slice from simulation T.OOL, T.IOL and T_40L. These slices are shown for a = 2, the plane wave begins to 
collapse at o = 1. Plane wave collapses along the vertical direction in these slices. The left panel is for T.OOL, middle panel is for T.IOL and the right panel 
is for T.40L. 



(T.lOP and T.40P) simulations there are no clumps but streams of 
particles that are falUng in and this reduces the number of scattering 
that take place - no « dependence means that dense streams run into 
each other head on with grazing collisions happening only rarely. 
Of course, in the simulation the presence of the plane wave leads to 
breaking of these streams into clumps as the streams are stretched 
inhomogeneously in the z direction. These clumps are aligned par- 
allel to the z axis. In the pancake region scattering of these streams 
occasionally leads to complex patterns. 

If presence of substructure and its growth in the pancake was 
the only cause for making the pancake thinner then pancake in these 
simulations should be thiimer as well. Figure 5 shows slices from 
simulations T_40L and T_40P for a = 2. A slice from the simu- 
lation PM.OOL is also plotted here for reference. This visual com- 
parison shows that the pancake is thinner in T_40L as compared to 
T.40P Indeed, the thickness of pancake in T.40P and PM.OOL is 
very similar. This reinforces the point that scattering of climips in 
the pancake region is the key reason for thinner pancakes. 

3.2 Pancakes and Viscosity 

The substructure is helping to confine the pancake to a smaller re- 
gion. It is interesting to study the collapse of a plane wave in an N- 
Body simulation and compare it with the collapse in the adhesion 
model (Gurbatov, Saichev and Shandarin, 1989) with a finite effec- 
tive viscosity. We first study the collapse of a plane wave in absence 
of any substructure, N-Body simulations PM.OOL for comparison 
with numerical simulations using the adhesion model with finite ef- 
fective viscosity. One dimensional adhesion simulations were done 
using the plane wave with the same amplitude as the N-Body simu- 
lations. We use the standard method for computing the trajectories 
of particles in the adhesion model (Weinberg and Gunn, 1990), a 
summary of the basic formalism is reproduced here for reference. 

In Adhesion approximation, the equation of motion for a par- 
ticle is replaced by the Burgers' equation (Gurbatov, Saichev and 
Shandarin, 1989; Weinberg and Gunn, 1990). In the one dimen- 
sional situation we are considering here, we have: 



du du 
db dx 



(5) 



Here u — dx/db is the 'velocity' of particles and b is the linear 
growth factor. This equation can be solved by introducing the ve- 
locity potential u = dip/dx, where ip coincides with the gravita- 
tional potential at the initial time. Solution has the following form. 
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dq. (7) 



Here q is the Lagrangian position of the particle and x is the Eule- 
rian position. In this method we integrate the differential equation 

for particle trajectories. At each time step velocity is calculated by 
above procedure at grid points and interpolated to particles posi- 
tions. 

Figure 7 shows the mass enclosed within a distance s from the 
centre of the pancake. The enclosed mass is defined as: 
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Here p{z) is the density at position z and Zc is the centre of the pan- 
cake. There is no ambiguity for comparing the results with N-Body 
simulations in case of no substructure as density depends only on z. 
While comparing other simulations with the adhesion solution, we 
will consider density averaged over x and y directions - Adhesion 
model is run only for the one dimensional problem. Top panel of 
figure 7 shows the enclosed mass M{z) at o = 2.0, middle panel 
is for a = 3 and the lower panel is for a = 4.0. The solid curve 
shows the enclosed mass for PM.OOL. In the region with a given 
number of streams, the N-Body curve is smooth. Jumps in mass 
enclosed occur at transition from single stream to multi stream re- 
gion, and at other transitions where the number of streams changes 
within the multi stream region. All other curves show A4{z) for 
adhesion model: dashed curve is for f = 400, dotted curve is for 
f = 600 and the dot-dashed curve is for v = 900. There is no 
constant effective viscosity curve that follows the N-Body curve 
closely through the multi stream regions. In regions with a given 
number of streams, the N-Body curve stays around a curve for con- 
stant effective viscosity in the adhesion model. A remarkable fact is 
that the N-Body curve for the three stream region at all the epochs 
follows the adhesion model curve for =~ 600. Similar behaviour 
is seen for the five stream region which follows v ~ 900 though 
the range of scales and epochs over which this can be resolved is 
somewhat limited. 

Addition of substructure clearly changes the character of the 
problem and the collapse is no longer one dimensional. However, 
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Figure 4. Top panel in this figure shows the density profile as a function of 
z, the direction of collapse for the plane wave. Density profile has been aver- 
aged over the directions transverse to the collapse of plane wave. The curves 
are for a = 2, simulations used are T.OOL (solid line), T_05L (dashed line), 
T.IOL (dot-dashed line), T.20L (dotted line) and T.40L (dot-dot-dashed 
line). The middle panel shows the rms transverse velocities of particles at 
the same epoch for T.05L (dashed line), T.IOL (dot-dashed line), TJ20L 
(dotted line) and T_40L (dot-dot-dashed line). The lower panel shows the 
rms transverse velocities of collapsed haloes at the same epoch for T_05L 
(dashed line), T.IOL (dot-dashed line), T.20L (dotted line) and T.40L (dot- 
dot-dashed line). 



the scale of the substructure is so small compared to the wavelength 
of the plane wave that the large scale collapse is still very close 
to planar. Figure 8 shows the mass enclosed within a distance s 
from the centre of the multi stream region for simulations PM_OOL, 
T.IOL and T_40L. Density is averaged over all x and y for this 
plot in the same manner as for figure 4. Also plotted in the figure 
are curves for the adhesion model (u = 600), where the calcula- 
tion is done without taking substructure into account. The motiva- 
tion for such a comparison is to see the effect of substructure on 
the favoured value of effective viscosity. Substructure removes the 
sharp change in density at the boundaries of 3-stream, 5-stream and 
7-stream regions and the curves for T.IOL and T_40L are smoother 
in the pancake region. The finite viscosity curve matches simu- 
lations with substructure over a wider range of scales than with 
PM.OOL. There are no other noteworthy differences. 



3.3 Mass Function 

Mass function of collapsed haloes in these simulations can be used 
to understand the influence of plane wave collapse on substructure. 
Collapsed structures form in these simulations primarily due to ini- 
tial density fluctuations at small scales, with some modulation by 
the collapse of the plane wave. In this section we study the effect of 
the collapsing plane wave on the mass function of collapsed haloes. 

These haloes were selected with the friends-of-friends (FOF) 
algorithm using a Unking length of I = 0.2 grid length. The initial 
power spectnmi has a peak at the scale corresponding to 1 /8 of the 
simulation box, or 16 grid lengths. Thus typical haloes will have a 
Lagrangian radius of about 8 grid lengths and should contain about 
500 particles. Thus a cutoff of 50 or more particles for haloes is 
reasonable for this study. 

In absence of the plane wave, the only perturbations are at 
small scales. The small scale perturbations are concentrated around 
a given mass scale and the mass function is also peaked around 
this mass at early epochs. At late epochs mergers lead to formation 
of haloes with a larger mass and the range of masses is greater 
for models with a larger amplitude of fluctuations. Figure 6 shows 
these features in the distribution of particles. These features can 
also be seen in figure 9 where mass fraction F{M) for a = 0.5, 
1.0 and 2.0 is plotted in different panels. F{M) is the fraction of 
total mass in collapsed haloes with halo mass above M . 

Adding the plane wave at a much larger scale than the small 
scale fluctuations essentially pushes much of the mass into the pan- 
cake region, leaving a small fraction of matter in the under dense re- 
gions that occupy much of the volume. Growth of small scale fluc- 
tuations in the under dense regions is inhibited whereas growth of 
fluctuations in the pancake region is enhanced, this is seen clearly in 
the slices from simulations shown in figure 6. Higher background 
density in the pancake region leads to rapid growth of perturba- 
tions, mergers of haloes also lead to formation of massive clumps. 
These effects become more pronounced at late epochs and result in 
a shift of mass function towards larger masses, indeed haloes at two 
distinct mass scales are present. Low mass clumps in under dense 
regions have the mass expected of haloes in regions where smaU 
scale fluctuation dominate whereas haloes of a much higher mass 
are present in the pancake region. Figure 9 shows these two mass 
scales very clearly. 

Total mass in collapsed haloes does not change significantly 
with the addition of the plane wave. Indeed for simulations T_40L 
and T_40, mass function is the same at a = 0.5 as small scales 
dominate. At late times (a = 2), the effect of plane wave makes the 
mass function of T.05L, T.IOL and T.20L similar. 

Not surprisingly, presence of large scale power leads to forma- 
tion of more massive haloes. However it does not seem to enhance 
the total mass in collapsed haloes. 



4 DISCUSSION 

In this paper we studied the effect of substructure on collapse of a 
plane wave. The key conclusions of the present study of the role of 
substructure are: 

• The pancake formed due to collapse of the plane wave is thin- 
ner if the in falling material is formed of collapsed substructure. 

• We show that collisions between clumps lead to enhancement 
of velocities transverse to the direction of large scale collapse. 

• We show that in simulations with substructure where colli- 
sions are suppressed, pancakes are not thinner. 
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Figure 5. This figure shows slices from simulations T_40L and T_40P for a = 2. The left panels shows a slice from the simulation PM.OOL, plotted here for 
reference. The central panel is for T_40P and the right panel is for T_40L. This visual comparison shows that the pancake is thinner in T_40L as compared to 
T_40P. Indeed, the thickness of pancake in T_40P and PM.OOL is very similar. 









Figure 6. This figure shows the effect of plane wave on evolution of small scale fluctuations. Panels in this figure show a slice from N-Body simulations. The 
top row is for a = 0.5, middle row is for a = 1 and the lower row is for a = 2. The left column is for T_05L, the second column is for T_05, the third column 
is for TJ20L and the right column is for T.20. 
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Figure 7. This figure sliows the mass enclosed within a distance z from the 
centre of the multi stream region. The top panel shows the curves for a = 2. 
The thick solid curve is for the N-Body simulation PM.OOL. Jumps in the 

mass enclosed occur at transition from multi stream region with 2ra + 1 
streams to 2n + 3 streams, with n a non-zero positive integer All other 
curves show M (z) for adhesion model: dashed curve is for u = 400, dotted 
curve is for f = 600 and the dot-dashed curve is for i/ = 900. The middle 
panel shows the same set curves for a = 3 and the lower panel is for a = 4. 



• Thus collision induced enhancement of motions transverse to 
the collapsing plane wave takes away kinetic energy from the di- 
rection of infall and leads to thinner pancakes. 

• Presence of large scale power shifts the mass function towards 
larger masses. There is, however, no change in the total mass in 
collapsed haloes. 

The points outlined above essentially relate to coupling of 
density fluctuations at well separated scales. Each of these points 
refers to a measurable effect of such a coupling. The nature of large 
scale fluctuation, a single plane wave, does not allow us to estimate 
the effect in terms of statistical indicators like the power spectrum. 
We plan to study these aspects with larger (256"^ ) simulations where 
the large scale collapse will also be generic. Large, high resolution 
studies are needed as 128^ simulations with particle mesh code 
have not shown any large effect in power spectrum at late times 
(Bagla and Padmanabhan, 1997a). 

Another important point to consider is that we have considered 
two well separated scales for fluctuations and there is no infall once 
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Figure 8. This figure shows the mass enclosed within a distance z from the 
centre of the multi stream region. The top panel shows the curves for a = 
2. The solid curve is for N-Body simulation PM.OOL. Other simulations 

are also plotted here T.IOL (dashed curve) and T_40L (dot-dashed curve). 
Dotted curve shows the mass enclosed in the one dimensional adhesion 
model with v = 600. The lower panel shows the same set of ciu-ves for 
a = 4 and the middle panel is for a = 3. 

fluctuations at the larger scales collapse. Numerical experiments 
that can shed light on effects of this feature are also required to 
improve our understanding of issues. 

We also compared the collapse of a plane wave in an N-Body 
with the collapse in the adhesion model with a finite effective vis- 
cosity. We found that: 

• The adhesion model predicts the variation of density very well 
with a constant effective viscosity in regions with a given number 
of streams. 

• Regions with a given number of streams coincide with the 
adhesion model with the same value of effective viscosity at all 
epochs. 



ACKNOWLEDGEMENTS 

JSB thanks R.Nityananda, T.Padmanabhan and K.Subramanian 
for useful discussions and suggestions. JSB also thanks Varun 
Sahni and Uriel Frisch for a useful discussion on related issues. 



Gravitational collapse and the role of substructure 9 




M 



Figure 9. This figure shows the cumulative mass function F{M) as a func- 
tion of mass M in N-Body simulations. The top panel is for a = 0.5, 
the middle panel is for a = 1 and the lower panel is for a = 2. Curves 
are shown for T_05 (solid curve), T_05L (thick solid curve), T_10 (dashed 
curve), T.IOL (thick dashed curve), T^O (dotted curve), T.20L (thick dot- 
ted curve), T_40 (dot-dashed curve) and T_40L (thick dot-dashed curve). 
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